source('lib.r')
#load('F3v2..RData')
# read conos object for T cells
scon=readRDS('Tcell.Normal.PTumor_new_conos.rds')
anoT = readRDS('final.anoT.rds')
scon$plotGraph(groups=anoT,size=0.1,plot.na=F,alpha=0.1)
set.seed(36)
anoT = as.factor(anoT)
anoT.pal <- setNames(rainbow(length(levels(anoT))),levels(anoT));
anoT.palf <- function(n) return(anoT.pal)
a2=scon$plotGraph(groups=anoT,raster=TRUE,plot.na=F,size=0.1,alpha=0.1,font.size = c(5, 5.5),palette=anoT.palf)
a2
#ggsave('F3a.new.pdf',a2,height=4,width=4)
# key marker gene expression
gs=c('CD3D','CD4','CD8A','FOXP3','GZMB','NKG7','FGFBP2','IL2', 'TNF' ,'IFNG', 'IL12A',"IL17F", "IL17A","RORC" ,'SELL','CCR7' )
lis=lapply(sn(gs),function(x) scon$plotGraph(gene=x,plot.na=F,size=0.1,alpha=0.3,title=x))
b= cowplot::plot_grid(plotlist=lis, ncol=4, nrow=4)
b
# dot plot for marker gene expression
library(cowplot)
anoT2 = anoT
gs=c('MKI67','TOP2A','CD8A','PDCD1','HAVCR2','GZMK','IFNG','GZMB','XCL1','XCL2','KLRG1','FGFBP2','NKG7','CMC1','NCAM1 ','KLRD1','NCR1','SELL','CCR7', "IL17A","RORC","CD4" , 'FOXP3', 'CTLA4','IL2RA')
#p2T=readRDS('Tcell.Normal.PTumor_new.p2Object.rds')
anoT2 = as.factor(anoT2)
anoT2=ordered(as.factor(anoT2),levels=c('Proliferation T','CTL-1','CTL-2','NKT','NK-1','NK-2','CD4 naive','Thelper','Treg'))
cname=names(anoT2)
aexp=t(p2T$counts)
cname = intersect(cname,colnames(aexp))
gs=intersect(gs,rownames(aexp))
#cname=intersect(cname,colnames(aexp))
p=Dotfig(gs,aexp[,cname],anoT2[cname],cols = c("blue","white", "red"))
p
# Exhaustion signature expression
gs=c('PDCD1','HAVCR2','FOXP3','CTLA4')
lis=lapply(sn(gs),function(x) scon$plotGraph(gene=x,plot.na=F,size=0.1,alpha=0.3,title=x))
b= cowplot::plot_grid(plotlist=lis, ncol=2, nrow=2)
b
# caculate exhausiton signature score
ylab='Exhaustion Score'
sp2= p2T
gs=readRDS('ExhaustionScore.csv')
cname=names(anoT)
cname=intersect(cname,rownames(sp2$counts))
df=Signature_score(anoT[cname],gs,sp2$counts,stype[cname],ssamp,min.num.cell=5,magnitude.normal = TRUE)
library(ggpubr)
df2=df[df$cell %in% c('CTL-2','CTL-1'),]
df2$cell=df2$fraction
tmp=df2
ysize=11
xsize=11
limHeight=1.2
name2=ylab
tmp$cell2=apply(data.frame(rownames(tmp)),1,function(x) strsplit(x,'[|]')[[1]][1])
p1 <- ggplot(tmp, aes(x=cell2,fill=fraction,y=score))+theme_classic() + geom_boxplot(outlier.shape = -1,width=0.5,position=position_dodge(width=0.7)) +theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + xlab("") + ylab(name2)
p1=p1+ geom_point(data = tmp,size=0.5,color=adjustcolor(1,alpha=0.3), position = position_jitterdodge(0.3))
p1=p1+ylim(c(min(tmp$score),max(tmp$score)*limHeight))
#p1=p1+ geom_point(data = tmp,color=adjustcolor(1,alpha=0.3),fill='grey', size = 1, shape = 21)
p1=p1+ theme(legend.position="none")
p1=p1+theme(axis.text.x = element_text(angle = 45, hjust = 1))
p1=p1+theme( axis.text.y = element_text(angle = 45, hjust = 0.5,color = "black"),axis.text.x=element_text(size=xsize,color = "black"),axis.title.y = element_text(size = ysize,color = "black"))
p1=p1+scale_fill_manual(values=fraction.palette1)
p1=p1+ scale_y_continuous(expand = c(0.04, 0.04), limits=c(min(df2$score), max(df2$score)*1.05))+ ggpubr::stat_compare_means(label = "p.signif", label.x = 1.5)
p1
# Treg activity
glist = scProcess:::getMarkers()
ylab='Treg activity'
gs=glist$TregActivity
cname=names(anoT)
cname=intersect(cname,rownames(sp2$counts))
df=Signature_score(anoT[cname],gs,sp2$counts,stype[cname],ssamp,min.num.cell=5,magnitude.normal = TRUE)
library(ggpubr)
#df2=df[df$fraction=='T-HG',]
drawBoxplot(ylab,df,ylab,myeloid.col=NULL,limHeight=1.2,height=2.7,width=2.4)
df2=df[df$cell=='Treg',]
df2$cell=df2$fraction
drawBoxplot('Treg',df2,ylab,myeloid.col=fraction.palette1,limHeight=1.08,height=2.4,width=2,dsize=2)
sig=compare_means(score ~ cell, data = df2)
sig
# cell proportions
cname=names(anoT)
ano2=data.frame('Cell'=anoT[cname],'SampleType'=ssamp[cname])
# Annotation vs sample
tmp2 <- acast(ano2, Cell ~ SampleType, fun.aggregate=length)
tmp3 <- (sweep(tmp2, 2, colSums(tmp2), FUN='/'))
tmp4 <- melt(tmp3)
head(tmp4)
names(tmp4) <- c('cell', 'sample','pc.of.sample')
tmp4$Group=NULL
tmp4$Group=sample.groups[as.character(tmp4$sample)]
p <- ggplot(na.omit(tmp4),aes(x=cell,y=pc.of.sample,dodge=Group,fill=Group))+geom_boxplot(notch=FALSE,outlier.shape=NA) + geom_point(position = position_jitterdodge(jitter.width=0.1),color=adjustcolor(1,alpha=0.3),pch=19,size=0.5)+theme_classic()+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_text(angle = 90, hjust = 0.5)) +xlab("") +ylab("fraction of total T cells")+theme(legend.position="top")+
scale_fill_manual(values=fraction.palette1)
p
library(ggpubr)
df=tmp4
rsig=NULL
for (i in unique(df[,1])){
tmp=df[df[,1]==i,]
sig=compare_means(pc.of.sample ~ Group, data = tmp) #
sig$cell=i
rsig=rbind(rsig,sig)
#sig[sig$p.signif!='ns',]
}
NKc = readRDS('/d0/home/meisl/Workplace/RCC/Figures/F3/conos/NK_conos.rds')
emb = readRDS('/d0/home/meisl/Workplace/RCC/Figures/F3/conos/NK.umap2.rds')
NKc$embedding = emb
gs = c('XCL1','CD44','FGFBP2','CX3CR1')
pl = lapply(sn(gs),function(x) NKc$plotGraph(gene=x,title=x,size=0.1,alpha=0.3,plot.na=F))
b= cowplot::plot_grid(plotlist=pl, ncol=2, nrow=2)
b
a2=NKc$plotGraph(groups=anof2,raster=TRUE,plot.na=F,size=0.4,alpha=0.4,font.size = c(5, 5.5),palette=anoT.palf)+theme(plot.margin = margin(0.1, 0.1, 0.1, 0.1, "cm"))
a2
# correlation of Proliferation T abundance and CTL-1 abundance
sample.groups
n1 = sample.groups[sample.groups=='Tumor'] %>% names()
n1
dat1=data.frame('x'=tmp3['Proliferation T',],'y'=tmp3['CTL-1',],'fraction'=sample.groups[colnames(tmp3)] )
dat1 = dat1[n1,]
ct <- cor.test(dat1$x,dat1$y,alternative='greater')
df=dat1
df$sample=rownames(df)
gg <- ggplot(na.omit(df), aes(x=x,y=y)) + geom_point(aes(color=fraction)) +
geom_smooth(method='lm',linetype='dashed',alpha=0.15,size=0.5) + theme_bw() +
theme(axis.text.y=element_text(angle=90)) + #guides(color=F) + #geom_text_repel(aes(color=sample))
xlab("Proliferation T abundance") + ylab("CTL-1 abundance") +scale_color_manual(values=fraction.palette1)+
geom_text(x=Inf,y=Inf,label=paste('R=',round(ct$estimate,2),' ','p=',round(ct$p.value,2),sep=''),hjust=1.2,vjust=1.2,size=3.5)
gg=gg+theme(legend.position='none')
gg
source('/d0-mendel/home/meisl/bin/FunctionLib/Lib/ploty.r')
#source('~/bin/violin.R')
modify_vlnplot.pvalue=function(gene,anoM,exp2,colp=NULL,pt.size=0){
cname=names(anoM)
ano2=anoM
#exp2=t(p2all$counts)
#data2=data.frame(exp2[gene,cname])
#colnames(data2)=gene
data2=data.frame(as.numeric(exp2[gene,cname]))
rownames(data2)=cname
colnames(data2)=gene
library(dplyr)
library(magrittr)
library(cowplot)
#a=SingleExIPlot.boxplot(
a=SingleExIPlot2(
data = data2,
idents = ano2,
split = NULL,
type = 'violin',
sort = FALSE,
y.max = NULL,
adjust = 1,
pt.size = pt.size,
cols =colp, # cols = myeloid.pal,
log = FALSE
)
#plot.margin = unit(c(-0.01, 0, -0.01, 0), "cm")
a=a+xlab("") + ylab(gene) + ggtitle("") +
theme(legend.position = "none",
# axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = rel(1), angle = 0),
axis.text.y = element_text(size = rel(1))#, plot.margin = plot.margin
)
df =SingleExIPlot2.dat(
data = data2,
idents = ano2,
split = NULL,
type = 'violin',
sort = FALSE,
y.max = NULL,
adjust = 1,
pt.size = pt.size,
cols =NULL, # cols = myeloid.pal,
log = FALSE
)
colnames(df)[1]='gene'
return(list('a'=a,'df'=df))
}
# HAVCR2 expressio in Tregs
cname = names(anoT[anoT=='Treg'])
p=modify_vlnplot.pvalue(gene='HAVCR2',as.factor(stype[cname]),t(p2T$counts),colp=fraction.palette1,pt.size=0.1)
a = p$a+theme( axis.title.y = element_text(angle = 90, hjust = 0.5,color = "black"))+scale_fill_manual(values = fraction.palette1)+xlab('')+theme(legend.position = 'none')+ theme(plot.title = element_text(hjust = 0.5))+ ggpubr::stat_compare_means(label.x = 1.5,label = "p.signif") #label = "p.signif",
df = p$df
a
a =a+ scale_y_continuous( expand=c(0, max(df$gene) * 0.1), limits=c(0, (max(df$gene) + max(df$gene) * 0.1)))